1 



Exact Local Bosonic Algorithm for Dynamical Quarks 

A.Galli a and Ph. de Forcrand 6 

a Max- Planck- Institut fur Physik, D-80805 Munich, Germany 
b SCSC-ETH-Zentrum, CH-8092 Zurich, Switzerland 



NO 
ON 
ON 



We present an exact local bosonic algorithm for the simulation of dynamical fermions in lattice QCD. We show 
that this algorithm is a viable alternative to the Hybrid Monte Carlo algorithm. 
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1. Introduction 

The local bosonic algorithm was proposed by 
Liischcr in a hermitian version as a new 
method to simulate full QCD. Here we present 
several significant improvements of this algo- 
rithm: an inexpensive stochastic Metropolis ac- 
cept/reject test to make the algorithm ex- 
act; a non-hermitian polynomial approximation 
(2|l| and a simple even-odd preconditioning . 
We illustrate the effectiveness of our algorithm 
for two-flavor QCD and compare its performance 
with the Hybrid Monte Carlo (HMC) algorithm 
in the relevant regime of small quark mass. 
We successfully model its static properties (the 
Monte Carlo acceptance), and try to disentangle 
its dynamics. From our analysis it appears that 
this version of the local bosonic algorithm is in- 
deed a viable alternative to the HMC. 

2. Description of the algorithm 

The full QCD partition function with two 
fermion flavors is given by 
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where D represents the fermion matrix and Sg 
the gauge action. This partition function can 
be approximated by introducing a polynomial ap- 
proximation of the inverse of the fermion matrix. 
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The polynomial P(z) = Yi 



Zk) of degree 
n is defined in the complex plane and approx- 
imates the inverse of z. Since we are investi- 



gating full QCD with two flavors, the determi- 
nant | det P(D)\ 2 manifestly factorizes into posi- 
tive pairs, so that the | dct prpyp term of the ap- 
proximation (|2|) can be expressed as a Gaussian 
integral over a set of boson fields 4>u {k = 1, n) 
with color and Dirac indices. The full QCD par- 
tition function ([!]) is then approximated by 
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where represents the set of all boson field fam- 
ilies, Sl is the local action Sl — Sg + Sb and 
S b = YJk=i \( D - z k)(/>k\ 2 - Making use of the lo- 
cality of Sl we may now simulate the partition 
function (^) by locally updating the boson fields 
and the gauge fields, using heat-bath and over- 
relaxation algorithms. 

The simulation of full QCD can be obtained from 
(|^) by correcting the errors due to the approx- 
imation through a Metropolis test at the end of 
each trajectory. Introducing the error term in the 
partition function we obtain 
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The correction term \dct(DP(D))\ 2 can 
be evaluated by expressing the determinant 
| det(DP(D))\ 2 as a Gaussian integral over an 
auxiliary field r\ and incorporating it in the par- 
tition sum (W) 



{d\J\\dn\\dtf\\d$\\d$\ 



by defining a new exact action 

S exact 
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where Sc(U,r]) — \[DP(D)]^ 1 rj\ 2 is the correction 
action. In carrying out the simulation one gener- 
ates configurations of (U, (f>) and rj such that the 
probability of finding a particular configuration 
is proportional to exp(— S eX act)- The strategy is 
to alternatively update the (U, (f>) fields and the 
77 field. The Metropolis acceptance probability 
P acc is defined using the exact action S exac t so 
that the transition probability of the algorithm 
satisfies detailed balance with respect to S exac t 
||. To be more precise: 

Paco{{UA)^{U',cl> l ))= (7) 

min (l, e-l^xl'+lxl 2 ^) if s G (U) > S G (U') 
min (l, e+l^'xlMxl 2 ^ if s G (U) < S G (U') 

where ^ is a Gaussian distributed spinor (used 
to generate rj with a global heat-bath by 
7] = DP(D)x or 77 = D'P(D')x) and W = 
[D'P(D')]- 1 DP(D). Here D and D' denote the 
fermion matrix with the U and U' gauge config- 
urations, 

respectively. The algorithm can be summarized 
as follows: 

• Generate a new Gaussian spinor \ with vari- 
ance one; 

• Update locally the boson and gauge fields m 
times (in a reversible and optimized order [Q) ac- 
cording to the approximate partition function (|^) ; 

• Accept/reject the new configuration according 
to the Metropolis acceptance probability (0). 

In order to evaluate the Metropolis acceptance 
probability we have to solve a linear system in- 
volving DP(D) or D'P(D'), for which we use the 
BiCGstab algorithm ||. This linear system is 
very well conditioned because P(D') (or P(D)) 
approximates the inverse of D' (or D). The cost 
for solving it is minimal and scales like the lo- 
cal updating algorithms in the volume and quark 
mass. 

3. Numerical results 

Numerical simulations using the exact local 
bosonic algorithm have been performed for dif- 
ferent lattice parameters. The majority of the 
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Figure 1 . Acceptance as a'runction of the number of 
fields. The line represents the fit ansatz (8). 

simulations are reported in ||. We explored the 
acceptance of the Metropolis correction test and 
the number of iterations of the BiCGStab al- 
gorithm used in that test for inverting DP(D). 
The study was performed by varying the degree 
n of the polynomial and the hopping parameter 
k. One observes that the acceptance increases 
quite rapidly with the degree of the polynomial, 
above some threshold || (see Fig.l). On the other 
hand, the number of iterations needed to invert 
DP(D) remains very low for high enough accep- 
tance. The data show clearly that the overhead 
due to the Metropolis test remains negligible pro- 
vided that the degree of the polynomial is tuned 
to have sufficient acceptance. Results obtained 
from simulations using even-odd preconditioning 
confirm, as expected, that the improvement of 
the approximation reduces the required number 
of bosonic fields by at least a factor two HQ]. 
This improvement is trivially implemented by a 
simple redefinition of the roots zj, ||] . 
We have also performed a simulation for param- 
eters relevant for physical interesting conditions. 
On a 8 4 lattice at (3 — 5.6 and k = 0.1585 (very 
near k c ) with Schrodinger functional boundary 
conditions we have compared the performance of 
our a 

lgorithm with HMC. In Table 1 we present the 
plaquette and its integrated autocorrelation time 
in units of fermion-matrix multiplications ([Dtp]) 
The simulation with the local bosonic algorithm 
was performed using 48 boson fields. Its average 
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acceptance is < P acc >= 0.589(1). The HMC 
data is obtained from These results show that 
for light quarks the peformance of our algorithm 
is comparable to the performance of HMC. 





< □ > 


Tint(O) 


HMC 
multiboson 


0.5726(1) 
0.57275(15) 


33750(6000) [D<f>] 
23040(4000) [D<t>] 



Table 1 



4. Optimisation 



Using some general assumptions and the known 
error bounds on the Chebyshev-like approxima- 
tion of -D -1 , we can obtain || an ansatz for the 
acceptance probability, which takes the form 
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where V is the lattice volume, k is the hopping 
parameter of the Wilson fermion matrix D and 
k c is the critical hopping parameter at a given 
(3. We expect the fitting parameter / to depend 
smoothly on f3, but very little on n, V and k. In 
|H we show that our ansatz (||) is very satisfac- 
tory (see Fig.l). At other values of /?, one can fix 
/ by a test on a small lattice, and then predict the 
acceptance for larger volumes and different quark 
masses. We can use our ansatz to analyze and op- 
timize the cost of a simulation. From the analysis 
it turns out that the optimal average acceptance 
is around 0.7-0.8. 

5. Scaling 

The total complexity with the volume and the 
quark mass is related to the scaling of the number 
of bosonic fields and the autocorrelation of the ob- 
servables. As the volume V increases, the number 
n of bosonic fields should grow like log V. As the 
quark mass m q decreases, the number of fields 
n must grow like m~ . This is a consequence 
of the exponential convergence of the polynomial 
approximation. We expect the autocorrelation to 
behave like r ~ nm~ ° IHg] so that the com- 
plexity is proportional to n 2 m~ a . The exponent 
a can only be derived from MC data since there is 
no clear theoretical understanding of the coupled 
dynamics of the gauge and boson fields. From an 



exploratory simulation j| on a 4 4 lattice at j3 = 
we obtained that a is near 1. This is only indica- 
tive, because the dependence of a on the volume 
and the coupling has yet to be explored. 

6. Conclusion 

We have presented an alternative algorithm to 
HMC for simulating dynamical quarks in lattice 
QCD. This algorithm is based on a local bosonic 
action. A non-hermitian polynomial approxima- 
tion of the inverse of the quark matrix is used 
to define the local bosonic action. The addition 
of a global Metropolis test corrects the system- 
atic errors. The overhead of the correction test is 
minimal. This algorithm is exact for any choice of 
the polynomial approximation. No critical tuning 
of the approximation parameters is needed. Only 
the efficiency of the algorithm, which can be mon- 
itored, will be affected by the choice of parame- 
ters. The cost of the algorithm scales similarly 
with that of HMC with the volume V of the lat- 
tice and with the inverse of the quark mass. Fi- 
nally, for light quarks the peformance of our algo- 
rithm is comparable to the performance of HMC 
provided that there is enough memory to store 
the needed boson fields. 
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